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ABSTRACT 



This thesis estimates the frequency response of a network where the only data is the 
output obtained from an Autoregressive-moving average (ARMA) model driven by a 
random input. 

Models of random processes and existing methods for solving AR.MA models are 
examined. The estimation is performed iteratively by using the Yule-Walker Equations 
in three different methods for the AR part and the Cholesky factorization for the MA 
part. The AR parameters are estimated initially, then MA parameters are estimated 
assuming that the AR parameters have been compensated for. After the estimation of 
each parameter set, the original time series is filtered via the inverse of the last estimate 
of the transfer function of an AR model or MA model, allowing better and better esti- 
mation of each model's coefficients. The iteration refers to the procedure of removing 
the .MA or AR part from the random process in an alternating fashion allowing the 
creation of an almost pure AR or MA process, respectively. As the iteration continues 
the estimates are improving. When the iteration reaches a point where the coefficients 
converge the last MA and AR model coefficients are retained as final estimates. 
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I. INTRODUCTION 



In this thesis, an iterative approach is presented to estimate the frequency response 
of a network where the only data is the output obtained from an Autoregressive-moving 
average (ARMA) model driven by a random input. The AR parameters are estimated 
initially, then the MA parameters are estimated assuming that the AR parameters have 
been compensated for.. 

To find the frequency response of an ARMA network two problems have to be ad- 
dressed. One is due to the shortness of the observed data (time limitation) and the re- 
sulting distortion of the estimated correlation function. The second problem is due to 
the nonlinear combination of the coefficients of the MA part as they appear in the esti- 
mate of the correlation function. 

To achieve this, we estimate coefficients of MA part and AR part iteratively. We 
use the Yule- Walker Method for AR coefficients estimation and the Cholesky decom- 
position for MA part estimation. We assume that the ARMA network output, i.e., the 
observed signal is produced by white Gaussian noise driving the network. Three differ- 
ent methods are used to estimate the AR coefficients. To minimize the effect of the 
MA part in the correlation function, correlation lags greater than the correlation length 
of the MA part are used in estimating the initial AR coefficients. The iteration refers 
to the procedure, which removes estimated .MA or estimated AR contribution from the 
output of the ARMA model in an alternating fashion. This allows better and better es- 
timates of the AR and MA coefficients as the iteration continues. For the iterative AR 
coefficient estimation three methods can be used, denoted by method 1, method 2 and 
method 3. .Method 1, 2 and 3 use correlation lags starting at p 4- 1, 0, and 0 respectively, 
where the first two methods use a square matrix inverse and the third method uses a 
Pseudo matrix inverse. 

Introduction to models of random processes is presented in Chapter II. Existing 
techniques to solve for the parameters of AR.MA models are presented in Chapter II. 
Chapter IV includes simulation results. Simulation studies employed the Matlab pack- 
age on the IB.M PC/AT. Conclusions and recommendations are presented in Chapter 
V. 
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II. MODELS OF RANDOM PROCESSES 



A. AR PROCESS MODEL 

Wc say that x(n) is an autoregressive process of order p, or simply an AR(p) process 
if it satisfies the dilference equation, 

,y(/i) + - 1 J + .... + -p) = t{n) (2.1) 

where dj, ..... are constant cocllicients. and e(;i) is a pure random process [Ref 1]. 
Equation (2.1) may be written in the following form, 

r 

A") = - -k) + e(/i) (2.2) 

A=1 



A realization of Eq.(2.2) is illustrated in Eigure 1. 



AR process x(n) 




^ 




(5 (J) 

-1 -1 


«- " ^ 





Figure 1. Autoregressive Model of Order p 



The system transfer function Il(z) between the input r.{n) and the output x(n) for 
the AR model shown in Figure 1 is 



2 



(2.3) 




It is required that A(z) has all its roots within the unit circle in the z-plane which 
guarantees that H(z) is a stable and causal filter. 

The form of equation (2.3) illustrates that AR models have finite poles but no ze- 
roes. Hence, this model is sometimes called an All-Pole model. Because an AR model 
has no poles outside or on the unit circle, it has the strict minimum-delay property and 
hence is always invertible. In general, minimum delay means that the transfer function 
must have no poles outside the unit circle, but can have poles on the unit circle. Strictly 
minimum delay means that the transfer function has no poles outside or on the unit 
circle. 

The AR model is also called an Infinite Impulse Response (IIR) filter. According 
to definition (2.2), output x(n) depends on past values of the output and on the present 
input. Because of this, it is also referred to as a Pure Feedback system. 

The power spectral density for AR models is given by 



and T is the sampling interval [Ref 2]. 

B. MA PROCESS MODEL 

The sequence x(n) is said to be a moving average process of order q (denoted by 
MA(q)) if it satisfies the difference equation. 




(2.4) 



where 



p 



/((fl = I + 



x{n) = bQt.[n) -i- b^z(n — 1) -h .... + b^e{n - q) 



(2.5) 



where b^, b ^, ...., ^,are coefficients, and e(«) is a pure random process [Ref 1). 
Equivalently, we may write. 



3 



( 2 . 6 ) 



9 

4") = X V(« - 

fc =0 



We may say that the output of an MA model depends only on present and past 
values of the input, i.e., there is no feedback in an MA model. 

A realization of Eq. (2.6) is illustrated in Figure 2. 




Figure 2. Moving Average Model of Order q 



The system transfer function H(z) between the input c{n) and the output x(n) for 
the MA process is 

7/(z) = D{z) = + b22~^ + •••• + (2.7) 

This transfer function has q finite zeroes, but no poles. Hence, the MA model is 
also called an All-Zero model. MA models are invertible if and only if B(z) has no zeroes 
outside the unit circle, nor on the unit circle. 

MA models are also called Finiie Impulse Response (FR) filters. 

The power spectral density for an MA model is given by 

= ( 2 . 8 ) 



where 



4 



<i 

£(/) = *0 + 

k=\ 

and T is the sampling interval. 

C. ARMA PROCESS MODEL 

We say that x(n) is an autoregressive-moving average process of order (p,q) or 
simply an ARMA(p,q) process if it is satisfies the difference equation, 

x(n) -f - 1) -f .... + cipX{n — p) = b^tin) + b^tin - 1) -f + b^t{n — p) (2.9) 

where, again, Uj, ...., a^, b ^, ...., b^ are coefficients and £(«) is a pure random process [Ref 
11 - 

Equation (2.9) may be written as. 



P oo 

- ^) + 

k=\ k=0 k=0 

The assumption b^ = 1 can be made without any loss of generality because the input 
e(«) can always be scaled to account for any filter gain [Ref 2]. 

A realization of Eq. (2.10) is illustrated in Figure 3. 

The system transfer function for the ARM.-\ process is given by 



^b^z{n - k) = Y^hi,c{n - k) 



( 2 . 10 ) 



H{z) = 



B{z) 

A(z) 



I + b^z ’ -t- £>2^ ^ -I- .... + b^z ^ 
1 -I- a,z~* + Q 2 Z~^ -t- .... + ctpZ~^ 



( 2 . 11 ) 



where A(z) is the z-transform of the AR part and B(z) is the z-transform of the MA part. 

Both polynomials A(z) and B(z) are assumed to have all of their zeros within the 
unit circle of the z-plane to guarantee that H(z) is a stable minimum-phase invertible 
filter. 

The power spectral density for ARMA models is given by [Ref 2] 



Barmaid — ^ 



1 

t 



B{f) 2 
AiJ) 



( 2 . 12 ) 
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D. RELATIONSHIPS OF RANDOM PROCESSES 

The Wold decomposition theorem (Ref. 3] relates the AR, MA and ARMA models. 
It shows that, if the Power Spectral Density is purely continuous, any AR or ARMA 
process can be represented by a unique MA model of infinite order. 

Another important theorem which is stated by Kolmogorov (Ref 4] says that any 
ARMA or MA process can be represented by an AR process of infinite order. 
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To illustrate these theorems, we model [Ref. 5] an ARMA(1,1) process by an 
AR(oo) or by an A/A(oo) process. From equation (2.1 1), the system transfer function for 
the ARM A( 1,1) process is 



//(z) = 



1 + ^)Z 
1 + 



If we use AR(oo) process to represent ARMA(1,1) process, where 

^(2) = 

1 + CjZ + C 2 Z + . . . 



and 



00 



\ 1 . \ -k ^ 

C(z) = 1 + / Ctz = 

^ 1+V 

k=] ‘ 

By using synthetic division we find that 

C(z) = 1 + (fl, — + {b\ — ai^,)z~^ + {b\ — a^b\)z~^ + 

Hence inverse z-transform of az"" is ad{k — m) [Ref 6] and 

'1 if A = 0 



d{k)^ 



0 else, 



the inverse z-transform of C(z) is 

Cfc = b{k) -F («, - b{)5{k -l) + (bt - a^b^)d{k -2) + {b\ - ay^)5{k - 3) + 



or 



_ j 1 if ^ = 0 



If we use a finite order AR(p) we should choose p to satisfy c^+i~0 or, equivalently 

Therefore, a high-order AR model will be required when the zero of the ARMA 
process gets closer to the unit circle. 
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In a similar way if we use an MA{oo) process to represent an ARMA(1,1) process, 



let 



H{z) = d(j cl^z * + d 2 Z ^ . 



where 




By using synthetic division as w'e did above, the inverse z-transform of D(z) will be 



If we use a finite order MA(q) model, we should choose q to satisfy or, 

equivalently 

Therefore, a high-order MA model will be required when the pole of the ARMA 
process gets closer to the unit circle. 

E. RELATIONSHIP OF AR, MA, AND ARMA PARAMETERS TO THE 
AUTOCORRELATION SEQUENCE 

In this section, we will present the relationship of the model parameters to the 
autocorrelation sequence [Ref 2]. 

If we multiply Eq. (2.10) by — m) and take expectation, the result will be 




1 if A: = 0 

(bi — a,)( — if A: > 1 




£{.Y(rt).v'X« - m)} = - - k)x*{n - m)) -b Y^b^E{t{n - k)x-{n - m)] (2. 13) 




where the superscript * is used to denote the complex conjugation. 
Equation (2.13) may be written as. 




(2.14) 



The cross correlation between the input and the output can be written as. 
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(2.15) 



e(« + 0 






/C=1 



^ ex (0 = '■„(0 + 

k=l 



where /zo = 1 by definition (Eq. 2.10). 

If we assume that the driving sequence is a white noise process of zero mean and 
variance erf then [Ref 2] 





0 


for 


/>0 




II 


2 


for 


/ = 0 


(2.16) 






for 


/<0 





When we substitute Eq. (2.16) in Eq. (2.14), we get final relationship between the 
ARMA parameters and the autocorrelation sequence. 



fxxM = 



kxx^i ~>n) for m<0 

p <? 

Yj^kPxxim -k) + olYj f'or 0<m<q 

p k=m 

- Y^(^k>'xx{>^ -k) for m > ? 



jfc=i 



(2.17) 



k=\ 



The relationship between the autocorrelation sequence and a pure autoregressive 
model may be written by setting ^ = 0 in Eq. (2.17) 



^xxi^) = 



- ^^k>'xxi’^ - k) for m > 



k=\ 



- -k) + a] for m = 0 

for m < 0 



(2.18) 



k=\ 

^xA -fn) 



The relationship between the autocorrelation sequence and a pure moving average 
model may be written by setting /> = 0 in Eq. (2.17) 
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0 



for m> q 









k=m 



>'xx*i -^) 



for 0 < m < q 
for m < 0 



In this case we should note that 

^k — ^k ^ ‘^k<q 



( 2 . 19 ) 
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III. TECHNIQUES TO SOLVE FOR THE PARAMETERS OF ARMA 

MODELS 



The estimation of the parameters of ARMA processes is a classical problem which 
is still being investigated by statisticians. Methods to find the parameters for purely AR 
processes are well known, but for ARMA processes some problems remain. 

There is a nonlinear relationship between the ARMA parameters and the 
autocorrelation of process x(n). The nonlinear equation (2.17) presents the difficulty of 
estimating the ARMA parameters, even when we know the autocorrelation sequence 
exactly. Techniques based on iterative maximum likelihood estimation (MLE) can be 
used to find the ARM.<\ parameters. These techniques require complex computations 
and are not guaranteed to converge, or they may converge to the wrong solution. 
Therefore they are not practical for real time series. For AR parameters, techniques 
based on the least squares criterion lead to solutions of linear equations and hence re- 
duce the computational complexity. Unfortunately, the moving average parameters of 
an ARMA model cannot be found easily by solving a set of linear equations. The MA 
parameters are convolved with the impulse response coefficients h(k) which causes a 
nonlinear relationship between the autocorrelation sequence and the filter coefficients. 

In section Ill-A and 111-B, we will discuss the methods of AR and M.'^ parameter 
estimation, and in section 111-C we will present an iterative approach to find the pa- 
rameters. 

A. AR PARAMETER ESTIMATION 

In this section, we present three widely used methods of extracting the model pa- 
rameters from a given block of measured data x(n). 

These methods are: 

1. The autocorrelation, or Yule-Walker method 

2. The covariance method 

3. Burg's method 

All three methods of estimating AR parameters are based on least-squares minimization 
criteria obtained by replacing the ensemble averages by appropriate time averages [Ref 

7 ]. 
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The criteria for the optimal forward ( e;(«) ) and backward ( e-{ri) ) predictors are 
obtained by minimizing 

£[c>)'] and £[cp»'] 

where ep{n) and e-{n) are the result of filtering x„ through the prediction-error filter as 
given by 



(^) T "I" ^Pi^n—l T + ^pp^n—p 

^p (^) ~ ^n-p T <^p^^n—p+\ T T ^pp^n 

The autocorrelation method is the most obvious and straightforward one. Equation 
(2.18) may be evaluated for the p + 1 lag indices 0^m<p and put into the following 
matrix form 



''xx(O) 


''xx(-l) 


rxxi -2) 


1 

1 


1 




2 

a. 


''jcx(l) 


'';c;c(0) 




''xx(-£+ 1) 


«i 




t 

0 






''xJc(O) 


-p + 2) 




= 


0 


>'xx{p) 


X 

. , 

1 


t-xxiP - 2) 


'•xx(O) 






0 



By noting that 

^xxi ~ ^xxi^) 

Eq. (3.1) can be written as 



^;c(0) 


''xx(l) 


rxxi2) 


rxx(p) 


1 




2 


''xx(l) 




''xx(l) 


>'xxiP - 1) 


a, 




t 

0 


''..(2) 


'■;cc(l) 




•'xxiP - 2) 


«2 


= 


0 


^xxip) 


•'xxip- 1) 


^xxip - 2) 


''xx(O) 


/P_ 







From Eq. (3.2) 

p 

= '^>'xx{‘)ai where ag = 1 

i=0 
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We can replace the autocorrelations r„(A) by the corresponding sample 
autocorrelations (biased estimate) computed from the given block of data, where 

X-\-k 

^ ^l+k^n Tor 0<k<p (3.3) 

«=0 



We use the biased autocorrelation estimator because the unbiased autocorrelation 
estimates may result in autocorrelation matrices that are not positive semidefinite, which 
means certain matrix equations have no solution. On the other hand, the autocorrelation 
matrices formed from the biased autocorrelation estimate will always be positive semi- 
definite [Ref 2 ]. 

By solving the normal equations, we can obtain estimates of the model's parameters 
{^ 1 , flj, flj,...., orf). Equation (3.1) is known as the AR, Yule-Walker or Normal 
Equations . The autocorrelation matrix of this equation is both Toeplitz and Hermitian 
because 



fx.x( ^xxi^) 



The solution of the Hermitian Toeplitz equations can be computed with the 
Levinson Algorithm. This algorithm is a lattice realization for linear prediction filters. 
We illustrate the prediction-error sequence in Figure 4. 

The prediction error is given by 

+ UpX„_^ + x„_2 + + appX„_p (3.4) 

Looking at the first two prediction errors 

e^{n) = x„->r + « 22-^«-2 

where (1, a,,), and (1, Oj,, a^^) represent the best predictors of orders p= 1 and p= 2, re- 
spectively. The extra index is used to indicate the order of the predictor. 

In the autocorrelation method the ensemble average in minimization of the forward 
prediction error is replaced by the least-squares time average criteria (p as follows 
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n=0 

A'+P-i (3.5) 

^J_ V 4-^ ^ 

) ■'^n + / ,‘^plc^n-k 

^ k=\ 

n=0 

where it is assumed that the data x„, .r,, , are observed and ,r„ = 0 for outside the 

range 0<n< A' — 1. The minimization of the time-average criteria with respect to the 

real and imaginary part of the 's will lead exactly to the same set of Yule-Walker 

Equations (3.1). One way to solve these equations is via the Levinson recursion which 

is an iterative technique that extracts the next order predictor from the previous one. 

1 he Levinson algorithm can be sunmiarized as follows (Ref. 7j: 

1- Initialize the recursion at p = 0, by setting 

Ao( 2 ) =1 and £q = '3cx(0) = 

where Ao(^) is the prediction-error filter and E„ is the mean-squared predic- 
tion error at the zero stage. So, initially we have no prediction. At stage p, 
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2 - 



the prediction-error filter will be Af,(z) and the mean-squared prediction er- 
ror will be E^. 

Compute which is called the reflection or PARCOR coefficient: 



p 




3- 



Recursively determine the {p -1- 1)'^ order prediction-error filter polynomial 






( 3 . 7 ) 



4- 



Update the mean-squared prediction error: 

£p*i - (I - 

Continue the iteration until the final desired order is reached. 



( 3 . 8 ) 



If the process x(n) is AR(p), then iteration will continue up to order p. It will pro- 
vide the AR coefficients a,,, a,, which are also the best prediction coefiicients. 

If we continue iteration after p, all prediction coefficients of order higher than p will be 
close to zero [Ref 7]. 

Although the autocorrelation method is the most obvious and efficient one, and the 
resulting prediction-error filter is guaranteed to be minimal phase, it suffers from the 
effect of prewindowing the data sequence x(n) by padding it with zeros to the left and 
to the right. This reduces the accuracy of the method, especially when we have short 
data records. 

In the Covariance method the ensemble average in minimization of the forward pre- 
diction error is replaced by the time average as follows 



The only difference between this method and autocorrelation method is in the limits 
of the summation. In the covariance method w'e observe all the data points needed to 
compute <p. Since we do not need the data outside the range 0<n<N— 1 we do not 




( 3 . 9 ) 
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need to set explicitely to zero. For that reason the covariance method seems to be 
more realistic. 

If we perform the minimization, we find the AR parameter estimates as the solution 
of equations (3.10) [Ref 5]. 

Cxx{p^) CxxM 

where 

Cxxii^k) = 

n=p 

Note that c^Jij,k) is an estimate of r„(j — k) . The c„(j,k) uses the sum of only N-p 
lag products to estimate the autocorrelation function for each lag even more lags are 
available. In contrast in the estimation of r„(0) the autocorrelation method uses all data 
point, while the covariance method uses only N-p data points in the summation. The 
minimization gives us a non-Toeplitz c„(j,k) matrix. This implies that we can not use 
the Levinson algorithm to solve Eq. (3.10). The equations may be solved by using the 
Cholesky decomposition which will be computationally more expensive. The estimated 
poles using this method are not guaranteed to lie within the unit circle. As an example, 
consider a first order predictor (p= 1) and a length-three (N= 3) sequence [Ref 7]. The 
time-average criteria will be 

2 

^ = y ^ I ej^(n) I ^ = y ((x, -I- Uj iXq)^ -I- (^2 + «i 

n=\ 

When we minimize $ by differentiating it with respect to a,, and setting the deriv- 
ative to zero , we have 



Cxxi^^P) 


■«r 




Cx;c(1.0) 


Cxxi^’P) 


«2 


= _ 


Cx;c(2,0) 


CxxM 






CxxiPfi) 



(3.10) 



A'-l 






(Xj -I- + (-^2 + «1 



a 



II ~ 



X^Xq -I- X2X1 

2 , 2 

Xq +X^ 



16 



Note the denominator does not depend on the variable jcj and if jtj is large enough 
we might have magnitude of a„ greater than one. In this case we do not get a 
minimum-phase prediction-error filter which is sometimes not desirable in practice. 

Burs's method estimates the reflection coefficients and then uses the Levinson 
recursion to estimate the AR parameters. Burg's minimization criteria is to minimize 
both the forward prediction error e‘{n) and the backward prediction error e~{n): 



A 

(p 



A-l 






n=p 



The computational steps are summarized below [Ref 7]: 



1 - 






Initialize, by setting ; 

e^{n) = eo"(n) = for 0 < « < A - 1 (3.11) 



and 



,v-i 




n=0 

At stage p- 1 , the prediction-error filter will be which is the Z trans- 
form of the sequence {L , . The mean-squared error 

will be £,_[. and e~_^{n) can be calculated for p — 1 < « ^ A' — 1. 



Compute the reflection coefficient: 



:V-1 

- 1 ) 

n=p 

r, = — (3.12) 



To guarantee the minimum-phase property, Y, must have a magnitude less 
than one. 



3- Compute A^{z) using the Levinson recursion given by: 



1 




1 




0 


^p,\ 








^p-\p-\ 


°p,1 


= 






^p-hp-2 






^p-'^,p-\ 










0 




1 



(3.13) 
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4 - 



Compute e^(n) and e~(n) for p < n < N — 1. 

^p(fj) = ^p-\(n) (3.14) 

^p(^) = Vi(" “ 1) “ 

5- Update the mean-squared error as follows: 

£,-(l-Y/£,_, (3.15) 

6- Continue the iteration until p equals the model order. 

The Burg's method estimates the poles which are on or inside the unit circle. This 
is due to the property | Y,| <1. Therefore, care must be taken to deal with the situation 
when I Y, I = 1, as this causes the prediction filter to become non-minimum phase. 

B. MA PARAMETER ESTIMATION. 

The most obvious approach to estimate the (MA parameters would be to solve the 
nonlinear equation (2.19) using the autocorrelation sequerice. Solutions of Eq. (2.19) 
involve difficult spectral factorization techniques [Ref 8]. 

There is another approach called Durbin s method which is related Kolmogorov's 
theorem [Ref 4] and is based on a high order AR approximation of the MA process. 
The AR process allows results using only linear operations. Let 

<t 

5(z)=l+^V"' 

*=i 

represent the system transfer function of an MA(q) process, and 



^00(2) = 1 + * 

*=i 

represent the system transfer function of an ARipo) process that is equivalent to the 
-MA(q) process. Therefore, we have 



B{2) = 



^00(2) 



or 



B{z)AJ,2)=\ 



The inverse Z-transform of the Eq. (3.16) is: 



(3.16) 
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(3.17) 






+ Yj^n 

n=\ 



= 5{m) = 



for 

for 



m = 0 

m = \,2,...,q 



where = 1 and = 0 for A: < 0. Therefore, the MA parameters can be determined from 
the infinite- order AR model by solving (3.17). 

In practice, one can calculate high-order AR(M) parameters, where M> q . Based 

on these parameter estimates (1, flv^(l), 4w(2), , an error in the MA part is 

computed [Ref 2]. 



9 

n=l 

According to Eq. (3.17) the error should be zero for all m except for m= 0. But in 
practice, the error will not be zero when using finite data, so MA parameter estimates 
are obtained by the minimization of the squared error power, given by 



P<i 



I 



I I 

M 



2 



m=0 



(3.19) 



This estimation procedure is an approximate maximum likelihood estimate (MLE). 
The approximate MLE procedure using Durbin's method for MA parameters results in 
the following estimation [Ref 5]. 



b = 



-R-'r 

'aa 



(3.20) 



where 




M- |/-y| 

A/+1 Zj 

n=0 



for /,j=l,2, ,q 



M-l 

[^■<3 = TTijrT 

«=0 



E 



^n^n+i 



for 



'■= 1, 2, ,q 
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Again the Levinson algorithm may be used to solve Eq. (3.20) for the b parameters. 
The estimated zeros of B(z) will be inside the unit circle by the minimum-phase property 
of the autocorrelation method. 

In summary, Durbin's method first uses the data jto, jt,, to find a large order 

AR(M) model using the autocorrelation method. Then using these AR parameter esti- 

A A A 

mates (1, a,, Am) as the data, (/?,, ^ 2 . •••, t>^) is found. 

Another technique for estimating the .MA parameters is to use the MA spectral 
estimator [Ref 5] which is given by Eq. (3.21) based on sample correlation values. 

(3-21) 

k=-q 

Since theoretically Eq. (3.21) is equal to a] \ the MA parameters can be found 

by using the Spectral factorization theorem [Ref 9]. This theorem shows that any ra- 
tional power spectral density of a stationary’ signal x„ can be factored into a minimum- 
phase form 

= (3.22) 

C. ITERATIVE ARMA PARAMETER ESTIMATION 

Let us consider the modelling of AR.MA transfer function as given in Figure 5. 

The computational steps are summarized below' for the iteration assuming know- 
ledge of the model order: 

1- Form the sample correlation of the time series y(n) by using the biased 
autocorrelation estimator. 

2- Get the AR coefficients using the correlation lags > q (to minimize the MA 
influence). 

3- Inverse filter the original time series y(n) via the inverse of the AR filter (use 
the AR coefficients of step 2) to get x^in) 

4- Form the sample correlation of the time series Xi{n) by using the biased 
autocorrelation estimator. 

5- Get the .M.A coefficients using the sample correlation of JC 2 (/t) . 

6- Inverse filter the original time series y(n) via the inverse of the MA filter 
(use MA coefficients w’hich are estimated in step 5) to get Xx{n). We assume 
that the M.A coefficients are of minimum phase (all roots are inside the unit 
circle). 
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AR(p) 


x^(n) . 


MA(q) 






W 


= ^ 


► y(n) 


£(n) ^ 














MA(q) 




AR(p) 








► y(n) 



Figure 5. Modelling of ARMA transfer function. 



7- 

8 - 
9- 



10 - 



11 - 

12 - 



13 - 



Form the sample correlation of time series ;c,(n) by using the biased 
autocorrelation estimator. 

Get the AR coeflicients using all lags or lags> q depending on the method 
used. 



Inverse filter the original y(n) series using the AR coefficient estimates 
which are obtained in step 8, to get .V 2 (/?). 

Get the sample correlation of by using the biased autocorrelation es- 
timator. 



Compute the MA coefficients using the sample correlation of .t 2 (n) . 



Compute the error in the AR and the MA parts as given below 

p ? 

enorij - 1) = + ^{b\ - for j = 2,3,.. .,10(3.23) 



!=l 



(=1 



where superscript j is used to denote the number of iteration. The upper 
count of j is e.xperimentally chosen. Iteration continue up to 10 if the co- 
efficients do not converge. 



If error > then go to step 6, else exit the program using the last updates 
of the filter coefficients. If j> 10 terminate with an error message. Note ). 
is a small e.xperimentally chosen number. 
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The AR coefTicients can be obtained via a pseudoinverse from the modified Yule- 
Walker equations [Ref. 2]. Equation (2.17) may be rewritten for the p lag indices 
q-\-\<m<q-\-p, and put into matrix form 



’'xx(‘}+ 1) 


rxxi^ - 1) 
rxxi^l) 


>'xxi(l-P + 1) 
>'xxi^ -p+2) 


'«r 




f'xJdl + I) 
^xxiq + 2) 


rxxi^l + p + 1) 


f'xxi^ + P + '2') 


1 






^xxi^l + P) 



We can compute the autocorrelation sequence for lags q-p+ 1 to q + p, therefore 
the AR parameters are found as the solution of Eq. (3.24) where the MA parts influence 
is minimal. The autocorrelation matrix is of Toeplitz form. 

The Moore-Penrose pseudoinverse /I' of an m x n matrix A provides the minimum- 
norm least squares solution to Ax = b. The solution is given by 

X = A^b 



where x is a n x I vector that simultaneously minimizes the squared equation error, for 
a given m x 1 vector b [Ref 2]. If m = n and rank A is n (i.e., A is nonsingular), then the 
pseudoinverse becomes the square matrix inverse A- = A~K If m> n (i.e., more 
equations than unknowns) and rank A is n, then 

A" = {A^A)~^A^ and 

x = (A^A)~'A^b 



The superscript H is used to denote the Hermitian transpose operation. This is the 
least squares solution for a set of overdetermined equations. 

In step 2 of the iteration, we will use correlation lags greater than q while in step 
8 either all lags or lags greater than q are used depending on the method. 

For MA coefficients estimation the Cholesky decomposition is used. If the matrix 
A is square and Hermitian, then the usual triangular factorization takes on the special 
form 



a = rr” 



(3.25) 



22 



where R is a lower triangular matrix with nonzero real principal diagonal elements. This 
decomposition is called the Cholesky decomposition [Ref. 2]. 

For the .VIA part, the statistical autocorrelation is given by 

q- \k\ 

««(*) - E 

n=0 

Let B be the lower triangular Toeplitz matrix where the matrix elements are given 
in terms of the impulse response of the filter B(z) [Ref 7]: 

^ni ~ ^n—i 

and let the autocorrelation matrix of x„ be 

Rxxi'j') = Rxxi‘ -J) 



Then, the transpose matrix will have matrix elements 



and Eq. (3.26) can be written as 



q- hWl 

Rxxi‘J) = l^xxi‘-J) = 

n=0 



RxxibJ) = 



Thus, in matrix notation 

r„=o]bb" 

which is related to Eq. (3.25) with the assumption that 




(3.27) 



Therefore, an approximation of the MA parameters can be found from the Cholesky 
decomposition of the correlation matrix i?„. 
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IV. SIMULATION RESULTS 



In this chapter, computer simulation results are presented to show how three dif- 
ferent estimation methods work on various ARMA models. Two ARMA(2,2) models 
which difier in pole-zero locations, an ARMA(2,3) model and an ARM.A(3,4) model are 
used as test models. In addition two realizations for each ARMA process are utilized. 
Two hundred data points are used in computing the sample autocorrelation values. The 
three dilTerent AR estimation methods are explained below. 

Method 1 

The .A.R coefficients are obtained via a square matrix inverse using the modified 
Yule- Walker equations. Correlation lags greater than q are used to minimize the M.-\ 
part influence at the first calculation. For the remainder of the iteration the same cor- 
relation lags are used to minimize potential influence from the MA part. A Cholesky 
factorization is used to find the MA part coefficients. 

Method 2 

The .AR coefficients are obtained via a square matrix inverse using the modified 
Yule-Walker equations. Correlation lags greater than q are used to minimize the MA 
part influence at first calculation. For the remainder of the iteration the correlation lags 
starting from zero are used assuming the MA part contribution has effectively been re- 
moved. The Cholesky factorization is used to find the MA part coefficients. 

Method 3 

The AR coefficients are obtained via a pseudoinverse instead of a square matrix 
inverse using the correlation function starting at the zero lag after first calculation. This 
allows the use of all important correlation lags. The Cholesky factorization is used to 
find the MA part coefficients. 

In addition, observation noise is added to one of the ARMA models and the re- 
sulting noisy sequence is processed via the three different methods. The observation 
noise is independent of the driving noise. It is white Gaussian noise with zero mean and 
unit variance. The signal to noise ratio (SNR) is about 15 dB. 

The computer program computes the differential errors in the AR and MA parts 
as given in Eq.(3.23), and then compares it with a small experimentally established value 
Z (i.e., A = 0.0001). If the error is less than 2 the program is terminated. If the error is 
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larger than ?. the program is reentered at step 6. Also, ifthe coeflicicnts do not converge 
the program will stop after nine iterations. 

Comparisons of the true and estimated coefficient differences, of pole-zero lo- 
cations, of distances between the true and estimated pole-zero locations and of radial 
differences between the true and estimated pole-zero locations are presented to show 
how well the three methods work. Also the spectra of the ARM A models are plotted by 
using the true and estimated coeflicients. 

A. THE ARMA(2,2) MODEL-A 

The pole-zero locations for this model are illustrated in Figure 6. 




Figure 6. The ARMA(2,2) Model-A, Pole-zero Locations 

I. METHOD I 

a. Noise Realization I 

The coefficients converge after two iterations. Tables 1 and 2 present the 

results. 
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Table 1. ARMA(2.2) MODEL-A, METHOD 1, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.80 


-1.8412 


-0.0412 


«2 


0.85 


0.8972 


-1-0.0472 




1.00 


1.3053 


4 - 0.3053 


bx 


0.80 


0.6439 


-0.1561 




0.80 


0.0572 


-0.7428 



Table 2. ARMA(2,2) MODEL-A, METHOD 1, POLE-ZERO COMPAR- 

ISON. REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.9 -1- 0.2j 


0.9206 -f- 0.2230j 


0.0308 


0.0189 


P2 


0.9 - 0.2j 


0.9206 - 0.2230] 


0.0308 


0.0189 




-0.4 + 0.8j 


-0.3771 


0.8003 


1.1071 


^2 


-0.4 - O.Sj 


-0.1162 


0.8488 


1.1071 



b. Noise Realization 2 

Using a different noise realization, the coefficients still converge after two 
iterations. Tables 3 and 4 present the results. 



Table 3. ARMA(2,2) MODEL-A, METHOD 1, 

COEFFICIENTS COMPARISON, 
REALIZATION 2 



Coeff. 


True 


Estimated 


Difference 




-1.80 


-1.8712 


-0.0712 




0.85 


0.9167 


+ 0.0667 


b. 


1.00 


1.1083 


+ 0.1083 


bx 


0.80 


0.5071 


-0.2929 


b2 


0.80 


0.1502 


-0.6498 
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Table 4. ARMA(2,2) MODEL-A, METHOD 1, POLE-ZERO COMPAR- 

ISON, REALIZATION 2 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.9 + 0.2j 


0.9356 + 0.2033j 


0.0358 


0.0047 


Pi 


0.9 - 0.2j 


0.9356 - 0.2033] 


0.0358 


0.0047 


Zl 


-0.4 4* 0.8j 


-0.2288 + 0.2884j 


0.5394 


0.2070 


^2 


-0.4 - 0.8j 


-0.2288 - 0.2884j 


0.5394 


0.2070 



The Spectra using the true and the estimated network coeflQcients are 
plotted in Figure 7. 

2. METHOD 2 



a. Noise Realization 1 



results. 



The coefficients converge after two iterations. Tables 5 and 6 present the 



Table 5. ARMA(2,2) MODEL-A, METHOD 2, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.80 


-2.0065 


-0.2065 


a. 


0.85 


1.0544 


+ 0.2044 


K 


1.00 


1.3866 


+ 0.3866 




0.80 


0.6467 


-0.1533 


bi 


0.80 


-0.0178 


-0.8178 
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Table 6. ARMA(2,2) MODEL-A, METHOD 2, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.9 + 0.2j 


1.0032 + 0.2189j 


0.1049 


0.0038 


P2 


0.9 - 0.2] 


1.0032 - 0.2189] 


0.1049 


0.0038 




-0.4 + 0.8j 


-0.4924 


0.8053 


1.1071 


^2 


0 

1 

p 


0.0261 


0.9064 


2.0344 



b. Noise Realization 2 

Using a different noise realization, the coefficients still converge after two 
iterations. Tables 7 and 8 present the results. 



Table 7. ARMA(2,2) MODEL-A, METHOD 2, 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 


«i 


-1.80 


-1.9267 


-0.1267 


«2 


0.85 


0.9729 


+ 0.1229 


b. 


1.00 


1.0982 


+ 0.0982 


bi 


0.80 


0.4598 


-0.3402 


b2 


0.80 


0.0902 


-0.7098 



Table 8. ARMA(2,2) MODEL-A, METHOD 2, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.9 + 0.2j 


0.9634 + 0.2118] 


0.0644 


0.0022 


Pi 


0.9 - 0.2j 


0.9634 - 0.21 18j 


0.0644 


0.0022 




-0.4 + 0.8j 


-0.2094 + 0.1958] 


0.6335 


0.3553 


^2 


-0.4 - O.Sj 


-0.2094 - 0.1958j 


0.6335 


0.3553 
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The spectra using the true and the estimated network coefTicients are plot- 
ted in Figure 8. 

3. METHOD 3 



a. Noise Realization 1 



results. 



The coefTicients converge after two iterations. Tables 9 and 10 present the 



Table 9. ARMA(2,2) MODEL-A, METHOD 3, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.80 


-1.8375 


-0.0375 


a. 


0.85 


0.8948 


+ 0.0448 


bo 


1.00 


1.3062 


+ 0.3062 




0.80 


0.6482 


-0.1518 


b2 


0.80 


0.0657 


-0.7343 



Table 10. ARMA(2,2) MODEL-A, METHOD 3, POLE-ZERO COM- 

PARISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.9 + 0.2j 


0.9187 + 0.2253j 


0.0314 


0.0218 


Pi 


0.9 - 0.2] 


0.9187 - 0.2253] 


0.0314 


0.0218 


Z\ 


-0.4 + 0.8j 


-0.3542 


0.8013 


1.1071 


h 


-0.4 - 0.8j 


-0.1421 


0.8405 


1.1071 



b. Noise Realization 2 

For a different noise realization, the coefficients stiU converge after one it- 
eration. Tables 11 and 12 present the results. 
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Table 11. ARMA(2.2) MODEL-A, METHOD 

3, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2 



Coeff. 


True 


Estimated 


Difference 




-1.80 


-1.8545 


-0.0545 




0.85 


0.9043 


-1-0.0543 


bo 


1.00 


1.1144 


+ 0.1144 


b: 


0.80 


0.5272 


-0.2728 


b2 


0.80 


0.1808 


-0.6192 



Table 12. ARMA(2,2) MODEL-A, METHOD 3, POLE-ZERO COM- 

PARISON, REALIZATION 2 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.9 + 0.2j 


0.9272 + 0.21 llj 


0.0293 


0.0051 


P2 


0.9 - 0.2j 


0.9272 - 0.21 llj 


0.0293 


0.0051 




-0.4 + 0.8j 


-0.2365 + 0.3260j 


0.5014 


0.1639 


^2 


-0.4 - 0.8j 


-0.2365 - 0.3260j 


0.5014 


0.1639 



The spectra using the true and the estimated network coefTicients are plot- 
ted in Figure 9. 

In summary’, the estimated coefTicients converge after two iterations in all 
three methods. For both noise realizations, the third method gives the best result in 
terms of the coefficients and pole-zero locations. For the noise realization 2, the coeffi- 
cients converge after one iteration using method 3. 

In each method the AR part coefficients of the ARMA(2,2) model tend to 
be more accurate than the MA part coefficients. 

B. THE ARMA(2,2) MODEL-B 

The pole -zero locations for this model are illustrated in Figure 10. 
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ARMA(2,2) Model A. Method 1. Realization 1 





Figure 7. The Spectra of ARMA(2,2) Model-A, Method 1 
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Figure 8. The Spectra of ARMA(2,2) Model-A, Method 2 
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Figure 9. The Spectra of ARMA(2,2) Model-A, Method 3 
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-1 - 0.5 0 0.5 1 



Figure 10. The ARMA(2,2) Model-B, Pole-zero Locations 



1. METHOD 1 

a. Noise Realizotioit I 

The coefficients converge after tliree iterations. Tables 13 and 14 present 



the results. 



Table 13. ARMA(2,2) MODEL-B, METHOD 



1, COEFFICIENTS COMPAR- 
ISON, REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.20 


-1.2614 


-0.0614 




0.52 


0.4143 


-0. 1057 


h 


1.00 


1.0518 


4-0.0518 




0.50 


0.4045 


-0.0955 


b2 


0.3125 


-0.1317 


-0.4442 
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Table 14. ARMA(2.2) MODEL-B, METHOD 1, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.6 -1- 0.4j 


0.6307 -1- 0.1285j 


0.2732 


0.3870 


Pi 


0.6 - 0.4j 


0.6307 - 0.1285j 


0.2732 


0.3870 


2, 


-0.25 + 0.5j 


-0.5950 


0.6074 


1.1071 


^2 


-0.25 - 0.5j 


0.2104 


0.6796 


2.0344 



b. Noise Realization 2 

For a difierent noise realization, the coeflicients converge after two iter- 
ations. Tables 15 and 16 present the results. 



Table 15. ARMA(2.2) MODEL-B, METHOD 

1, COEFFICIENTS COMPAR- 
ISON. REALIZATION 2 



Coeff. 


True 


Estimated 


Difference 




-1.20 


-1.6704 


-0.4704 




0.52 


0.7879 


+ 0.2679 


K 


1.00 


1.1018 


+ 0.1018 


h, 


0.50 


0.1188 


-0.3812 


Ih 


0.3125 


-0.3508 


-0.6633 



Table 16. ARMA(2,2) MODEL-B, METHOD I, POLE-ZERO COMPAR- 

ISON, REALIZATION 2 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4] 


0.8352 + 0.3006] 


0.2553 


0.2425 


Pi 


0.6 - 0.4j 


0.8352 - 0.3006] 


0.2553 


0.2425 




-0.25 + 0.5] 


-0.6207 


0.6224 


1.1071 


-2 


-0.25 - 0.5j 


0.5129 


0.9121 


2.0344 
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The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 11. 

2. METHOD 2 

a. Noise Realization I 

The coefficients converge after six iterations. Tables 17 and 18 present the 

results. 

Table 17. ARMA(2,2) MODEL-B, METHOD 



2, COEFFICIENTS COMPAR- 
ISON, REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.20 


-1.1531 


-1-0.0469 


CI2 


0.52 


0.5013 


-0.01S7 


h 


1.00 


1.0174 


-H 0.0 174 


b. 


0.50 


0.4732 


-0.0268 


b: 


0.3125 


0.1202 


-0.1923 



Table 18. ARMA(2,2) MODEL-B, METHOD 2, POLE-ZERO COMPAR- 



ISON. REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


/’i 


0.6 + 0.4j 


0.5766 + 0.41 lOj 


0.0258 


0.0312 


P2 


0.6 - 0.4j 


0.5766 - 0.41 lOj 


0.0258 


0.0312 


^1 


-0.25 -t- 0.5j 


-0.2326 -F 0.2532] 


0.2474 


0.2793 


Z2 


1 

P 

C/i 

1 

P 


-0.2326 - 0.2532] 


0.2474 


0.2793 



b. Noise Realization 2 

For a different noise realization, the coefficients converge after five iter- 
ations. Tables 19 and 20 present the results. 
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Table 19. ARMA(2.2) MODEL-B, METHOD 

2, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2 



Coeff. 


True 


Estimated 


Difference 


a, 


-1.20 


-1.5914 


-0.3914 




0.52 


0.7540 


+ 0.2340 


bo 


1.00 


1.0610 


+ 0.0610 


b, 


0.50 


0.1301 


-0.3699 


b. 


0.3125 


-0.2894 


-0.6019 



Table 20. ARMA(2,2) MODEL-B, METHOD 2, POLE-ZERO COMPAR- 

ISON. REALIZATION 2 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4j 


0.7957 + 0.3477] 


0.2025 


0.1760 


P2 


0.6 - 0.4j 


0.7957 - 0.3477] 


0.2025 


0.1760 


Zj 


-0.25 + 0.5j 


-0.5872 


0.6030 


1.1071 


^2 


-0.25 - 0.5j 


0.4645 


0.8720 


2.0344 



The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 12. 

3. METHOD 3 



a. Noise Realization 1 

The coefficients converge after four iterations. Tables 21 and 22 present 

the results. 
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Table 21. ARMA(2.2) MODEL-B, METHOD 

3. COEFFICIENTS COMPAR- 
ISON, REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 


a, 


-1.20 


-1.1818 


-1-0.0182 


a. 


0.52 


0.5245 


-H 0.0045 


h 


1.00 


1.0030 


-1-0.0030 


b. 


0.50 


0.4436 


-0.0564 




0.3125 


0.1015 


-0.2110 



Table 22. ARMA(2,2) MODEL-B, METHOD 3, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 -1- 0.4j 


0.5909 + 0.4187] 


0.0207 


0.0284 


P: 


0.6 - 0.4j 


0.5909 - 0.4187] 


0.0207 


0.0284 


^1 


-0.25 4 - 0.5] 


-0.2211 + 0.2286] 


0.2729 


0.3050 


Z 2 


-0.25 - 0.5j 


-0.2211 - 0.2286] 


0.2729 


0.3050 



b. Noise Realization 2 

For a dilTerent noise realization, the coeflTicients converge after two iter- 
ations. Tables 23 and 24 present the results. 



Table 23. ARMA(2,2) MODEL-B, METHOD 

3, COEFFICIENTS COMPAR- 
ISON, REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 


a, 


-1.20 


-1.5105 


-0.3105 




0.52 


0.6962 


+ 0.1762 


h 


1.00 


1.0488 


+ 0.0488 


bi 


0.50 


0.1827 


-0.3173 


h 


0.3125 


-0.2154 


-0.5279 
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Table 24. ARMA(2.2) MODEL-B, METHOD 3, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4j 


0.7553 + 0.3547j 


0.1617 


0.1489 


P2 


0.6 - 0.4j 


0.7553 - 0.3547j 


0.1617 


0.1489 


h 


-0.25 + 0.5j 


-0.5486 


0.5823 


1.1071 


22 


-0.25 - 0.5j 


0.3744 


0.7999 


2.0344 



The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 13. 

In summary, under each noise realization the method 3 gives the best result 
in terms of the coefficients and pole-zero locations. 

For each method, the spectrum of the ARMA(2,2) model using estimated 
coefficients indicates reasonably accurate pole locations but the zero locations of the 
spectrum do not follow the original ones. This means that the AR part coefficients of 
the ARMA model tend to be more accurate than the MA part coefficients. 

C. THE ARMA(3,4) MODEL 

The pole-zero locations for this model are illustrated in Figure 14. 

1. METHOD 1 

a. Noise Realization 1 

The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 
in the comparison. Tables 25 and 26 present the results. 
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AUMA(2,2) Model U. Method 1, Henlizalion I 
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Figure 11. llie Spectra of ARMA(2,2) Model-B, Method 1 
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AUMA(2.2) Model U. Method 2. Healizotloii 1 





Figure 12. The Spectra of ARMA(2,2) Model-B, Method 2 
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Figiiie 13. The Spectra of ARMA(2,2) Model-B, Method 3 
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Figure 14. The ARMA(3,4) Model Pole-zero Locations 



Table 25. ARMA(3,4) MODEL, METHOD 1, 
COEFFICIENTS COMPARISON, 
REALIZATION I. 



Coelf. 


line 


Estimated 


Dilfcrence 


<T| 


-1.60 


-1.5275 


+ 0.0725 


<'2 


2.18 


1.5642 


-0.6158 


<h 


-1.36 


-1.0165 


+ 0.3435 


(h 


0.7225 


0.4656 


-0.2569 




1.00 


3.3059 


+ 2.3059 




0.40 


1.3473 


0.9473 


Ih 


0.48 


-1.6482 


-2.1282 


l-h 


-0.32 


-2.4187 


-2.0987 
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Table 26. ARMA(3,4) MODEL, METHOD 1, POLE-ZERO COMPAR- 

ISON. REALIZATION I. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.2 -(- 0.9] 


0.0718 + O.SlOlj 


0.1565 


0.1302 


Pi 


0.2 - 0.9j 


0.0718 - 0.8101] 


0.1565 


0.1302 


Pi 


0.6 -1- 0.7] 


0.6919 + 0.4745j 


0.2435 


0.2610 


Pi 


O 

• 

o 


0.6919 - 0.4745] 


0.2435 


0.2610 




-0.4 + O.Sj 


-0.6754 + 0.5652j 


0.3619 


0.4103 


^2 


O 

o 

1 


-0.6754 - 0.5652] 


0.3619 


0.4103 


^3 


0.4 


0.9433 


0.5433 


0.0000 



b. Noise Realization 2 

Using a different noise realization, the coefficients still do not converge in 
nine iterations. Because of oscillations about two values, the average of the two values 
is used as an estimate of the coefficients in the comparison. Tables 27 and 28 present 
the results. 



Table 27. ARMA(3,4) MODEL, METHOD 1, 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-2.0479 


-0.4479 


a. 


2.18 


2.2690 


+ 0.0890 


«3 


-1.36 


-1.4743 


-0.1143 


^4 


0.7225 


0.9468 


+ 0.2243 




1.00 


2.6176 


+ 1.6176 


b. 


0.40 


0.8737 


+ 0.4737 


bi 


0.48 


-1.2400 


-1.7200 


bi 


-0.32 


-1.7033 


-1.3833 
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Table 28. ARMA(3.4) MODEL, METHOD I, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Px 


0.2 + 0.9j 


0.0518 + 0.8258] 


0.1657 


0.1560 


P2 


0.2 - 0.9] 


0.0518 - 0.8258] 


0.1657 


0.1560 


Pi 


0.6 + 0.7] 


0.9722 + 0.6618] 


0.3741 


0.2644 


Pi 


0.6 - 0.7j 


0.9722 - 0.6618j 


0.3741 


0.2644 


h 


-0.4 + 0.8j 


-0.6316 + 0.5489] 


0.3415 


0.3916 


h 


-0.4 - 0.8j 


-0.6316 - 0.5489] 


0.3415 


0.3916 


h 


0.4 


0.9294 


0.5294 


0.0000 



The spectra using the true and the estimated network coefTicients are plot- 
ted in Figure 15. 



2. METHOD 2 



a. Noise Realization 1 



results. 



The coefTicients converge after nine iterations. Tables 29 and 30 present the 



Table 29. ARMA(3,4) MODEL, METHOD 2, 

COEFFICIENTS 

COMPARISON.REALIZATION I. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-1.6186 


-0.0186 




2.18 


2.2093 


+ 0.0293 




-1.36 


-1.3948 


-0.0348 


a. 


0.7225 


0.7366 


+ 0.0141 


t>o 


1.00 


1.1463 


+ 0.1463 


bx 


0.40 


0.2664 


-0. 1336 


b2 


0.48 


-0.0262 


-0.5062 


bi 


-0.32 


-0.3386 


-0.0186 
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Table 30. ARMA(3,4) MODEL, METHOD 2, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.2 -1- 0.9j 


0.2033 -f- 0.9077j 


0.0083 


0.0016 


P2 


0.2 - 0.9j 


0.2033 - 0.9077j 


0.0083 


0.0016 


Pi 


0.6 -1- 0.7j 


0.6060 -1- 0.6957j 


0.0073 


0.0079 


A 


0.6 - 0.7j 


0.6060 - 0.6957j 


0.0073 


0.0079 




-0.4 -f- O.Sj 


-0.4197 -f- 0.5572j 


0.2435 


0.1819 




1 

0 

1 

o 


-0.4197 - 0.5572j 


0.2435 


0.1819 


^3 


0.4 


0.6070 


0.2070 


0.0000 



b. Noise Realization 2 

Using the second noise realization, the coefficients converge after nine it- 
erations. Tables 31 and 32 present the results. 



Table 31, ARMA(3,4) MODEL, METHOD 2. 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-1.6668 


-0.0668 


^2 


2.18 


2.3645 


+ 0. 1 845 




-1.36 


-1.5282 


-0.1682 




0.7225 


0.8527 


-1-0.1302 


bo 


1.00 


0.9923 


-0.0077 


b, 


0.40 


0.1274 


-0.2726 


bi 


0.48 


0.2222 


-0.2578 


bi 


-0.32 


-0.1221 


-1-0.1979 
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Table 32. ARMA(3,4) MODEL, METHOD 2, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.2 + 0.9j 


0.2063 + 0.9315] 


0.0321 


0.0007 


P2 


0.2 - 0.9j 


0.2063 - 0.93 15j 


0.0321 


0.0007 


ft 


0.6 + 0.7j 


0.6271 0.7372j 


0.0460 


0.0037 


ft 


0.6 - 0.7j 


0.6271 - 0.7372] 


0.0460 


0.0037 




-0.4 + O.Sj 


-0.2286 + 0.5674] 


0.2889 


0.0806 


Z2 


-0.4 - O.Sj 


-0.2286 - 0.5674] 


0.2889 


0.0806 


^3 


0.4 


0.3287 


0.0713 


0.0000 



The spectra using the true and the estintated network coefTicients are plot- 
ted in Figure 1 6. 



3. METHOD 3 



a. Noise Realization I 



results. 



The coefTicients converge after six iterations. Tables 33 and 34 present the 



Table 33. ARMA(3,4) MODEL, METHOD 3, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-1.5821 


-0.0179 


C7. 


2.18 


2.1352 


-0.0448 


ft 


-1.36 


-1.3329 


+ 0.0271 


ft 


0.7225 


0.7023 


-0.0202 


bo 


1.00 


1.1733 


+ 0.1733 


bi 


0.40 


0.3140 


-0.0880 


h 


0.48 


-0.0604 


-0.5404 


b, 


1 

o 


-0.3978 


-0.0778 



47 




Table 34. ARMA(3,4) MODEL. METHOD 3, POLE-ZERO COMPAR- 



ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.2 -1- 0.9j 


0.1913 -1- 0.9j 


0.0087 


0.0092 


Pi 


0.2 - 0.9j 


0.1913 -0.9j 


0.0087 


0.0092 


Ps 


0.6 -1- 0.7j 


0.5997 -1- 0.6855j 


0.0145 


0.0101 


A 


0.6 - 0.7j 


0.5997 - 0.6855j 


0.0145 


0.0101 


2i 


-0.4 O.Sj 


-0.4533 -1- 0.5691] 


0.2369 


0.2089 


Z2 


-0.4 - O.Sj 


-0.4533 - 0.5691] 


0.2369 


0.2089 




0.4 


0.6406 


0.2406 


0.0000 



b. Noise Realization 2 

Using the second noise realization, the coefficients converge after four it- 
erations. Tables 35 and 36 present the results. 



Table 35. ARMA(3.4) MODEL, METHOD 3. 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-1.5671 


+ 0.0329 


^2 


2.18 


2.2062 


+ 0.0262 




-1.36 


-1.4005 


-0.0405 


(h 


0.7225 


0.7728 


+ 0.0503 


bo 


1.00 


0.9932 


-0.0068 




0.40 


0.2049 


-0.1951 


^2 


0.48 


0.2119 


-0.2681 


b, 


-0.32 


-0.1654 


+ 0.1546 
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Table 36. ARMA(3,4) MODEL, METHOD 3, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.2 + 0.9j 


0.1876 + 0.9336j 


0.0358 


0.0203 


Pi 


0.2 - 0.9j 


0.1876 - 0.9336j 


0.0358 


0.0203 


Pi 


0.6 H- 0.7] 


0.5959 + 0.7051] 


0.0065 


0.0069 


A 


0.6 - 0.7] 


0.5959 - 0.7051] 


0.0065 


0.0069 


Z] 


-0.4 4- O.Sj 


-0.2937 + 0.5923j 


0.2333 


0.0033 


Z: 


-0.4 - 0.8] 


-0.2937 - 0.5923] 


0.2333 


0.0033 


Z 3 


0.4 


0.3810 


0.0190 


0.0000 



The spectra using the true and the estimated network coefficients are plot- 
ted in Figure 17. 

When using method 1 the coefficients do not converge. They oscillate 
about two sets of values in an alternating fashion, hence the average of the two sets is 
used as the estimate of the coefficients in the comparison. The parameters of the first 
realization converge after nine iterations using method 2, and after six iterations using 
method 3. The parameters of the second realization converge after four iterations using 
method 3. Method 3 provides the best results in terms of the coefficients and pole-zero 
locations for both noise realizations. The coefficients converge also earlier when using 
method 3 compared to the other two methods. Method 2 performs not as well but gives 
also reasonable results. 

The spectrum due to the poles of the ARMA(3,4) models using the esti- 
mated coefficients closely resembles the original spectrum. The AR part coefficient es- 
timates are more accurate than the MA part coefficient estimates. 

D. THE ARMA(3,4) MODEL WITH OBSERVATION NOISE 

Observation noise is added to the ARM A( 3,4) model and resulting noisy sequence 
is processed via the three different methods. The observation noise is independent of the 
driving noise. It is white Gaussian noise with zero mean and unit variance. The signal 
to noise ratio is about 15 dB. 

This model is illustrated in Figure 18. 
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AUMAO.-I) Model. Method 1. Realization 1 





Figure 15. Hie Spectra of ARI\IA(3,4) Models, Method I 
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ARMA(3,4) Model, Method 2, Realization 1 





Figure 16. Tlie Spectra of ARMA(3,4) Models, Method 2 
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AUMA(3.4) Model, Melliod 3. Realization 1 





Figure 17. The Spectra of ARMA(3,4) Models, Method 3 
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Figure 18. Tlie ARMA(3,4) Mode! with Observation Noise 



1. METHOD 1 

The coeH'icients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coellicients 
in the comparison. Tables 37 and 38 present the results. 



Table 37. ARMA(3,4) MODEL WHH OB- 

SERVATION NOISE, METHOD 1, 
COEFFICIENTS COMPARISON. 



Coelf. 


True 


Estimated 


Dillerence 




-1.60 


-2.0285 


-0.4285 




2.18 


1.9241 


-0.2559 




-1.36 


-1.4285 


-0.0685 




0.7225 


0.8145 


H- 0.0920 




1. 00 


4.5732 


-t- 3.5732 


h 


0.40 


1.7315 


1.3315 




0.48 


-2.4570 


-2.9370 


l-h 


-0.32 


-3.3648 


-3.0448 
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Table 38. ARMA(3,4) MODEL WITH OBSERVATION NOISE, 

METHOD 1. POLE-ZERO COMPARISON 



Poles- 

Zeios 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.2 + 0.9j 


0.0147 + 0.8350j 


0.1963 


0.2010 


P2 


0.2 - 0.9j 


0.0147 - 0.8350j- 


0.1963 


0.2010 


Pi 


0.6 + 0.7] 


0.9995 + 0.4109] 


0.4931 


0.4721 


A 


P 

CS 

1 

p 


0.9995 - 0.4109] 


0.4931 


0.4721 




-0.4 + 0.8] 


-0.6723 + 0.5565j 


0.3652 


0.4157 


^2 


o 

p 


-0.6723 - 0.5565j 


0.3652 


0.4157 


^3 


0.4 


0.9660 


0.5660 


0.0000 



The spectra using the true and the estimated network coefficients are plotted 
in Figure 19. 

2. METHOD 2 

The coefficients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefficients 
in the comparison. Tables 39 and 40 present the results. 



Table 39. ARMA(3,4) MODEL WITH OB- 

SERVATION NOISE, METHOD 2, 
COEFFICIENTS COMPARISON. 



Coeff. 


True 


Estimated 


Difference 




-1.60 


-2.6027 


-1.0027 


CI 2 


2.18 


3.8865 


+ 1.7065 


A 


-1.36 


-2.8006 


-1.4406 


A 


0.7225 


1.4933 


-1-0.7708 


bo 


1.00 


8.5983 


+ 7.5983 


bi 


0.40 


3.4234 


+ 3.0234 


b2 


0.48 


-3.2534 


-3.7334 


bi 


-0.32 


-4.6219 


-4.3019 
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Table 40. ARMA(3,4) MODEL WITH OBSERVATION NOISE, 



METHOD 2, POLE-ZERO COMPARISON 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.2 + 0.9j 


0.3264 -1- 0.8586j 


0.1330 


0.1446 


P2 


0.2 - 0.9j 


0.3264 - 0.8586] 


0.1330 


0.1446 


Pi 


0.6 + 0.7] 


0.9749 + 0.9052] 


0.4273 


0.1138 


A 


0.6 - 0.7j 


0.9749 - 0.9052] 


0.4273 


0.1138 




-0.4 -1- 0.8j 


-0.6152 -t- 0.5170] 


0.3555 


0.4082 


^2 


-0.4 - O.Sj 


-0.6152 - 0.5170] 


0.3555 


0.4082 


^3 


0.4 


0.8323 


0.4323 


0.0000 



The spectra using the true and the estimated network coefficients are plotted 
in Figure 20. 

3. METHOD 3 

The coefficients converge after four iterations. Tables 41 and 42 present the 

results. 



Table 41. ARMA(3,4) MODEL WITH OB- 

SERVATION NOISE. METHOD 3, 
COEFFICIENTS COMPARISON. 



Coeff. 


True 


Estimated 


Difference 


«i 


-1.60 


-1.4994 


+ 0.1006 


«2 


2.18 


2.1311 


-0.0489 


«3 


-1.36 


-1.3224 


+ 0.0376 


^4 


0.7225 


0.7353 


+ 0.0128 




1.00 


1.3898 


+ 0.3898 


bi 


0.40 


0.2875 


-0.1125 


bi 


0.48 


0.0298 


-0.4502 


b. 


-0.32 


-0.2462 


+ 0.0738 
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Table 42. ARMA(3,4) MODEL WITH OBSERVATION NOISE, 



METHOD 3, POLE-ZERO COMPARISON 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.2 + 0.9j 


0.1823 -H 0.9331j 


0.0375 


0.0257 


Pi 


0.2 - 0.9] 


0.1823 - 0.9331] 


0.0375 


0.0257 


Pi 


0.6 + 0.7] 


0.5674 + 0.7010] 


0.0326 


0.0281 


Pa 


0.6 - 0.7j 


0.5674 - 0.7010]' 


0.0326 


0.0281 




-0.4 -t- 0.8] 


-0.3481 -H 0.4908] 


0.3135 


0.1532 


^2 


O 

• 

o 

1 


-0.3481 -0.4908] 


0.3135 


0.1532 


^3 


0.4 


0.4893 


0.0893 


0.0000 



The spectra using the true and the estimated network coefficients are plotted 
in Figure 21. 

In sununars’, the estimated coefficients do not converge using method 1 or 2 
but they do converge after four iterations using method 3. The spectra using method 1 
and 2 are poor. But the method 3 gives results close to the actual ones. The spectrum 
estimated with this method follows the original pattern except for the zero location. 
This is believed to be partially due to the imprecise MA coefficient estimation procedure. 
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Figure 19. llie Spedra of ARMA(3.4) Models uilh Obser. Noise, Method 1 
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AUMA(3,'I) Model willi Obser. Noise. Method 2 




Figure 20. The Spectra of ARMA(3,4) Models uith Obser. Noise, Method 2 
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A!U!A(3,i) Model with Obser. Noise, Method 3 




Figure 21. Tlie Spectra of ARMA(3,4) Models uitli Obser. Noise, Method 3 
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V. CONCLUSIONS AND RECOMMENDATIONS 



In this thesis, an iterative technique to estimate the coeflicients of ARMA models 
driven by a random input is presented. The AR parameters are estimated initially, then 
MA parameters are estimated assuming the AR parameters have been compensated for. 
In an iterative fashion MA and AR contributions are removed from the original data 
allowing improved AR and MA coefficient estimates. Three different AR estimation 
methods are experimentally explored. The third method provides the best result in terms 
of the true and estimated coefficients and pole-zero locations. This method uses the 
pseudoinverse with correlation function values starting at the zero lag after removing the 
VIA influence. For models of real time data which have an odd number of poles, one 
should address the issue of the DC component (i.e., remove it). Simulation results for 
an odd ordered pole model is presented in Appendix C. 

Also, by examining the spectra of the models, we can say estimation of poles is 
obtained more accurately than the estimation of zeros. The Cholesky factorization is 
used for the estimation of the .MA part coefficients but it is an approximate solution. 
For that reason, we do not expect superior results for the VIA estimation part. 

Further research should concentrate on improving the .MA part coefficients esti- 
mation. 



60 



APPENDIX A. COMPUTER PROGRAM TO ESTIMATE 

COEFFICIENTS 



^*^V**')Wc*‘>V?V:>V*Vc'5Wr5V*Vf'5V**‘5V'5Wc*')Wf ********* 

% This program written by * 

% Gurhan Kayahan for IBM PC/AT * 

% Matlab package. * 

^******Vf*********************Vc*****Vc* 

% GENERATION OF ARMA TIME SERIES * 

4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4\ 4% 4\ 4\ 4\ 4% 4\ 4% 4\ 4\ 4\ 4 C 4 \ 4\ 4\ 4\ 4\ 4% 4\ 

i=sqrt(-l); 

rootsl=[ -0. 25 + 0. 5i - 0. 25 - 0. 5i]; change zeros for each model 

roots2= [0, 6 4 - 0. 4i 0, 6 - 0. 4i 0. 8]; change poles for each model 

b=poly( roots 1); find true coefficients of MA part 

a=poly( roots2); f ind true coefficients of AR part 

a=real(a); 

t=l: 1000; 

y=t*0; 

rand( ^normal ^ ); generate random signal 
stddev=sqrt( 1); 
rnum=s tdev*rand( t ) ; 

yl=filter(b, a, rnum); filter random signal through ARMA model 
x= [ones(l , 999)]; generate step signal 

y2=filter(b,a,x); find step response of ARMA model to discard transient response 
y=yl(52:252); 

V'f Vv ^ V V? V? Vc V'T Vc Vr Vc Vc Vr Vc Vc Vc* V? Vr V? V? Vc V? Vc Vc V? V* Vc Vc Vr Vc V* V* Vr Vc* Vc 

% CALCULATE AR COEFFICIENTS USING YULE -WALKER APPROACH * 

rl=xcorr(y); 
rl=rl(201: 399); 
r=flipy(rl); 
row=r( 194; 196); 
col=flipy(r(192: 194)); 
a=toepl itz( col , row) ; 
b=flipy(-r(191: 193)); 
c=b'; 

x=b*inv(a); square matrix inverse 
al=x' ; 

V'T Vv Vr Vv V\ V'T VrVcVc V'T /V /V^V Vc Vc VcVrVc Vc V% Vc Vc « 

% CALCULATE MA COEFFICIENTS USING CHOLESKY FACTORIZATION * 



d= [1]; 

n= [1 al(l) al(2) al(3)]; 
z=filter(n,d,y); 
zl=z(4; 201); 
rzl=xcorr(zl); 
rz=rzl(198; 395); 
i=l; 

while i<199 

rz3( i)=rz( i)/198; 
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i=i+l; 

end 

rz2=rz3( 1: 3); 
rzt=toeplitz( rz2); 
cho=chol(rzt); 
bl=cho( 1 5 : ) ; 

VcVr *'}Wr’}Wr>V'iWc"}V /V* *** ■jV’jWoV'sV'jV'jV •>V*')V’}V*'}V*')V'>V'5V'}W«’‘3V'>Wc -sWc 



% ROUTINE FOR ESTIMATION OF AR COEFFICIENTS * 

^V?'sV'>VVc’}V’}V’}VVf’}W«’^V’}V')V’5V">V*’)V'jWr*^V')V"}Ws’TV'}V'}V"5V'5V"3VVoWoW«"jV'>Wr’}Hc’**’}V'A’'>V 

i=2; 

while i<ll 

d= 1,1) bl{i- 1,2) bl{i- 1,3)]; 

n= [1]; 

zar=filter(n,d,y); 
rar=xcorr ( zar ) ; 
rl=rar(201: 399); 
r=flipy(rl); 

row=r( 194: 196); change correlation lags for each method 
col=flipy(r( 192: 194)); change correlation lags for each method 
g=toeplitz(col , row); 

e=flipy( -r( 191: 193) ); change correlation lags for each model 
f=e^ 



x=f*inv(^); square matrix inverse (replace with x=pinv(g) for pseudo inverse) 

a(i, : )=x ; 

al= Car, a(i,: )']'; 

^‘V5V:V*Vr'}VVrVc’V?Vr*')Wc’Vr:VVf’}WcVc*'>Wr’}W?V«'^V'jWoV’jWoV*V?’;V:iWc”>Wc')W«'*’}V?V'>Wf 

% ROUTINE FOR ESTIMATION OF MA COEFFICIENTS 



^‘V*Vc*Vc’*VcVrVcVcVv***VrVf*V*Vr‘>‘V*Vc'>V‘V‘5V'>V*VoV*VfV>‘VvVv’***Vr^V*>V*4V^V'5V* 



d= [1]; 

n= [1 a(i,l) a{i,2) a(i,3)]; 

z=filter(n,d,y); 

zl=z(4: 201); 

rzl=xcorr(zl); 

rz=rzl( 198: 395); 

k=l; 

while k<199 



rz3(k)=rz(k)/198; 

k=k+l; 

end 

rz2=rz3( 1: 3); 
rzt=toeplitz( rz2); 
cho=chol( rzt); 
bt(i,: )=cho(l,: ); 

bl= [br, bt(i,: )']'; 

e(i-l)=sum( ( al( i , : ) -al( i- 1 , : ) ) . a 2 )/3+sum( (bl( i, : ) -bl( i-1 , : ) ). a 2 )/3; 
if e(i-l) < 0. 0001 
al=al; 
bl=bl; 

else 



i=i+l; 

end 

end 
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APPENDIX B. COMPUTER PROGRAM TO CALCULATE TRUE AND 

ESTIMATED SPECTRA 



X CALCULATION OF TRUE SPECTRUM 



a= [1 “1.20 0.52]; change true poles for each models 

b= [1 0.50 0.3125]; change true zeros for each models 

al= la zeros(l, 123)]; zero padding 

bl= lb zero.s(l , 124)]; zero padding 

fal=fft(al); fast fourier transform of AR part 

fbl=fft(bl); fast fourier transform of MA part 

mfal=abs( fal), a 2; 

mfbl=abs(fbl). a 2; 

c=mfbl. /mfal; 

a=max(c); 

i=l; 

while i<129 
c(i)=c(i)/a; 
i=i+l; 



end 

spec=10*logl0( c); 



% CALCULATION OF ESTIMATED SPECTRUM 



%**■ 



' Vf ‘V Vf Vf * ‘V * * * Vf !>V Vc Vc if * Vr i: Vc * ‘V !>V “V ■ 



■VvVf 






a= [1 “1.1531 0.5013]; change estimated poles for each models 

b= [1.0174 0.4732 0.1202]; change estimated zeros for each models 

al= la zeros(l , 123)]; 

bl= lb zeros(l , 124)]; 

fal=fft(al); 

fbl=fft(bl); 

mf al=abs( f al) . a 2; 

mfbl=abs(fbl). a 2; 

c=mfbl. /mfal; 

a=max(c); 

i=l; 

while i<129 
c(i)=c(i)/a; 
i=i+l; 



end 

spec7=10*logl0( c); 



O/ 

/©*' “ '' '' 



. y . y y . y . ^ ^ y y ^ ^ ^ 



% PLOTTING OF TRUE AND ESTIMATED SPECTRUM 



^ Vc V'T Vv V'T Vc 



t=0; 127; 
t=t'; 

plot(t(l: 65),spec7(l;65),'*' ,t(l; 65),spec(l: 65)) 
title( ' ARMA(2 ,2) Model, Method 3, Realization l') 
xlabel( ' Sample Number') 
ylabel( 'Magnitude(dB) ' ) 
grid 

text(25,-5,' Using true coefficients') 
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text( 25 , -9 , ' **** After 4 iterations ' ) 
delete sped, met 
meta sped 
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APPENDIX C. SIMULATION RESULTS FOR AN ARMA(2,3) MODEL 



For models with an odd number of poles, one of the poles must be located on the 
real axis in the z-domain to generate real time data. 

1 he pole-zero locations for this model are illustrated in Figure 22. 




Figure 22, The ARMA(2,3) Model Pole-zero Locations 

I. METHOD 1 

a. Noise Realization I 

The coefTicients do not converge in nine iterations. Because of oscillations 
about two values, the average of the two values is used as an estimate of the coefFicients 
in the comparison. Tables 43 and 44 present the results. 
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Table 43. ARMA(2,3) MODEL, METHOD 1, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-2.00 


-1.9425 


-t- 0.0575 




1.48 


1.4028 


-0.0772 




-0.416 


-0.3639 


-h 0.0521 


bo 


1.00 


1.0562 


-1-0.0562 


b, 


0.50 


0.5184 


-0.0184 


^2 


0.3125 


0.1681 


-0.1444 



Table 44. ARMA(2,3) MODEL, METHOD 1, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4j 


0.6526 + 0.3809j 


0.0559 


0.0596 


P2 


0.6 - 0.4j 


0.6526 - 0.3809j 


0.0559 


0.0596 


Pi 


0.8 


0.6373 


0.1627 


0.0000 




-0.25 + 0.5j 


-0.2454 + 0.3145j 


0.1855 


0.1989 




-0.25 - 0.5j 


-0.2454 - 0.3145j 


0.1855 


0.1989 



b. Noise Realization 2 

Using a difierent noise realization, the coefTicients converge after seven it- 
erations. Tables 45 and 46 present the results. 
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Table 45. ARMA(2,3) MODEL, METHOD 1, 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-2.00 


-2.8571 


-0.8571 




1.48 


2.7525 


1.2725 




-0.416 


-0.9162 


+ 1.3322 


bo 


1.00 


1.3418 


+ 0.3418 


b. 


0.50 


0.1291 


-0.3709 


b2 


0.3125 


-0.4401 


-0.7526 



Table 46. ARMA(2.3) MODEL, METHOD 1, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4j 


0.8301 + 0.2765J 


0.2611 


0.2664 


P2 


0.6 - 0.4j 


0.8301 - 0.2765j 


0.2611 


0.2664 


Pi 


0.8 


1.1968 


0.3968 


0.0000 




-0.25 + 0.5j 


-0.6229 


0.6237 


1.1071 




-0.25 - 0.5j 


-0.5266 


0.9236 


2.0344 



2. METHOD 2 

a. Noise Realization 1 

The coefTicients do not converge in nine iterations. Because of oscillations 
about two values, the average of the fw'o values is used as an estimate of the coefficients 
in the comparison. Tables 47 and 48 present the results. 
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Table 47. ARMA(2.3) MODEL, METHOD 2, 

COEFFICIENTS COMPARISON, 
REALIZATION I. 



Coeff. 


True 


Estimated 


Difference 


a, 


-2.00 


-1.6534 


+ 0.3466 




1.48 


0.9727 


-0.5073 




-0.416 


-0.2044 


+ 0.2116 


bo 


1.00 


1.2677 


+ 0.2677 


b, 


0.50 


0.8274 


+ 0.3274 


bi 


0.3125 


0.3741 


+ 0.0616 



Table 48. ARMA(2,3) MODEL, METHOD 2, POLE-ZERO COMPAR- 

ISON. REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.6 + 0.4j 


0.5267 + 0.2515j 


0.1656 


0.1425 


P2 


0.6 - 0.4j 


0.5267 - 0.25 15j 


0.1656 


0.1425 


Pi 


0.8 


0.6001 


0.1999 


0.0000 


*^1 


-0.25 + 0.5j 


-0.3263 + 0.4349] 


0.1006 


0.1807 


"2 


-0.25 - 0.5j 


-0.3263 - 0.4349j 


0.1006 


0.1807 



b. Noise Realization 2 

Using a second noise realization, the coeflicients do not converge. Because 
of oscillations about two values, the average of the two values is used as an estimate of 
the coefficients in the comparison. Tables 49 and 50 present the results. 
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Table 49. ARMA(2,3) MODEL, METHOD 2, 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-2.00 


-1.1005 


-t- 0.8995 




1.4S 


-0.1144 


-1.5944 




-0.416 


0.3945 


+ 0.8105 


bo 


1.00 


1.8387 


+ 0.8387 


b, 


0.50 


1.3064 


+ 0.8064 


bi 


0.3125 


0.4440 


4- 0.1 3 15 



Table 50. ARMA(2,3) MODEL, METHOD 2, POLE-ZERO COMPAR- 

ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


P\ 


0.6 + 0.4] 


0.S145 + 0.2SS3j 


0.2418 


0.2478 


P2 


0.6 - 0.4j 


0.8145 - 0.2SS3j 


0.2418 


0.2478 


/b 


0.8 


-0.5285 


1.3285 


3.1415 




-0.25 + 0.5] 


-0.3553 + 0.3395j 


0.1919 


0.3444 


^2 


-0.25 - 0.5] 


-0.3553 - 0.3395] 


0.1919 


0.3444 



3. METHOD 3 



a. 



results. 



Noise Realization 1 

The coefl'icients converge after five iterations. Tables 51 and 52 present the 



Table 51. ARMA(2,3) MODEL, METHOD 3, 

COEFFICIENTS COMPARISON, 
REALIZATION 1. 



Coeff. 


True 


Estimated 


Difference 




-2.00 


-1.1057 


-f- 0.8943 




1.48 


0.2111 


-1.2689 


«3 


-0.416 


0.0890 


+ 0.5050 


bo 


1.00 


1.8594 


+ 0.8594 


b, 


0.50 


1.4885 


-0.9885 


b. 


0.3125 


0.8320 


-0.5195 



Table 52. ARMA(2,3) MODEL, METHOD 3, POLE-ZERO COMPAR- 

ISON, REALIZATION 1. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


/’i 


0.6 + 0.4] 


0.6493 + 0.2] 


0.2059 


0.2892 


P: 


0.6 - 0.4] 


0.6493 - 0.2j 


0.2059 


0.2892 


P', 


0.8 


-0.1928 


0.9928 


3.1415 




-0.25 + 0.5] 


-0.4002 + 0.5360] 


0.0509 


0.1777 




-0.25 - 0.5j 


-0.4002 - 0.5360j 


0.0509 


0.1777 



b. Noise Realization 2 

Using a second noise realization, the coefficients do not converge. Because 
of oscillations about two values, the average of the two values is used as an estimate of 
the coefficients in the comparison. Tables 53 and 54 present the results. 
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Table 53. ARMA(2,3) MODEL, METHOD 3, 

COEFFICIENTS COMPARISON, 
REALIZATION 2. 



Coeff. 


True 


Estimated 


Difference 




-2.00 


-1.7973 


-1-0.2027 




1.48 


0.9959 


-0.4841 




-0.416 


-0.9354 


-0.5194 


b. 


1.00 


1.2540 


-1- 0.2540 


bi 


0.50 


0.6150 


-1-0.1150 


bi 


0.3125 


0.0151 


-0.2974 



Table 54. ARMA(2,3) MODEL, METHOD 3, POLE-ZERO COMPAR- 



ISON, REALIZATION 2. 



Poles- 

Zeros 


True 


Estimated 


Distance 


Radial Diff. 


Pi 


0.6 + 0.4j 


0.1263 + 0.7679j 


0.5997 


0.2564 


Pi 


o 

1 

\o 

o 


0.1263 - 0.7679] 


0.5997 


0.2564 


Pi 


0.8 


1.5446 


0.7446 


0.0000 




-0.25 + 0.5j 


-0.4645 


0.5440 


1.1071 


^2 


-0.25 - 0.5j 


-0.0259 


0.5479 


1.1071 



Simulation results have shown that for the ARMA(2,3) model with one 
pole on positive real axis, realization 2 gives poor results. This is because the DC com- 
ponent for this realization is - 0.9017 while it is -0.0580 for realization 1. For parameter 
estimations the DC component should be removed. Therefore the data is filtered effec- 
tively removing the pole on the positive real axis (i.e., DC component is removed and 
the ARMA(2,3) model becomes ARMA(2,2) model). The results of this ARMA(2.2) 
model was presented in chapter IV-B, The AR.MA(2,3) and AR.MA(2,2) true and esti- 
mated spectra using method 1 are presented in Figure 23. Method 2 results are pre- 
sented in Figure 24 and method 3 results are presented in Figure 25. The simulation 
shows that the method 3 gives good results as long as the DC component is small. 
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At(MA(2.n) MndrI. Mrltiud t. I 



AUMA(2.Z) tl. MrUiiid t. rtmtlrnllnM I 






Figure 23. The Spectra of ARMA(2,3) and ARMA(2,2) Model-B, Metliod I 
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Figure 24. 1 lie Spectra of ARMA(2,3) and ARMA(2,2) Model-B, Metliod 2 
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